#
# step8 (PCA) Perform linear dimensional reduction
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))

# Examine and visualize PCA results a few different ways
#print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
# VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
DimPlot(pbmc, reduction = "pca")

# DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
# DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)





()RunPCA
./seurat-4.1.0/R/generics.R:440:RunPCA <- function(object, ...) {


#' Run Principal Component Analysis
#'
#' Run a PCA dimensionality reduction. For details about stored PCA calculation
#' parameters, see \code{PrintPCAParams}.
#'
#' @param object An object
#' @param ... Arguments passed to other methods and IRLBA
#'
#' @return Returns Seurat object with the PCA calculation stored in the reductions slot
#'
#' @export
#'
#' @rdname RunPCA
#' @export RunPCA
#'
RunPCA <- function(object, ...) {
  UseMethod(generic = 'RunPCA', object = object)
}






()
./seurat-4.1.0/R/dimensional_reduction.R:976:RunPCA.Seurat <- function(


#' @param reduction.name dimensional reduction name,  pca by default
#'
#' @rdname RunPCA
#' @concept dimensional_reduction
#' @export
#' @method RunPCA Seurat
#'
RunPCA.Seurat <- function(
  object,
  assay = NULL,
  features = NULL,
  npcs = 50,
  rev.pca = FALSE,
  weight.by.var = TRUE,
  verbose = TRUE,
  ndims.print = 1:5,
  nfeatures.print = 30,
  reduction.name = "pca",
  reduction.key = "PC_",
  seed.use = 42, #默认随机数种子
  ...
) {
  #如果为空，则使用默认 assay
  assay <- assay %||% DefaultAssay(object = object)
  # 获取 assay 对象
  assay.data <- GetAssay(object = object, assay = assay)

  # 把任务转给 Assay 对象的同名方法
  reduction.data <- RunPCA(
    object = assay.data,
    assay = assay,
    features = features,
    npcs = npcs,
    rev.pca = rev.pca,
    weight.by.var = weight.by.var,
    verbose = verbose,
    ndims.print = ndims.print,
    nfeatures.print = nfeatures.print,
    reduction.key = reduction.key,
    seed.use = seed.use,
    ...
  )
  
  #> pbmc0@reductions #覆盖前
  #list()
  # 覆盖到 pca 这个slot中
  object[[reduction.name]] <- reduction.data
  #覆盖后 
  #> pbmc@reductions
  #$pca
  #A dimensional reduction object with key PC_ 
  # Number of dimensions: 50 
  # Projected dimensional reduction calculated:  FALSE 
  # Jackstraw run: FALSE 
  # Computed using assay: RNA

  # 记录日志
  object <- LogSeuratCommand(object = object)
  return(object)
}















()
./seurat-4.1.0/R/dimensional_reduction.R:933:RunPCA.Assay <- function(

没有指定 features 则使用高变基因。
features 必须在 scaled data 中，如果不在或者 variance=0 则舍弃，用其余的做PCA。

#' @param features Features to compute PCA on. If features=NULL, PCA will be run
#' using the variable features for the Assay. Note that the features must be present
#' in the scaled data. Any requested features that are not scaled or have 0 variance
#' will be dropped, and the PCA will be run using the remaining features.
#'
#' @rdname RunPCA
#' @concept dimensional_reduction
#' @export
#' @method RunPCA Assay
#'
RunPCA.Assay <- function(
  object,
  assay = NULL,
  features = NULL,
  npcs = 50,
  rev.pca = FALSE,
  weight.by.var = TRUE,
  verbose = TRUE,
  ndims.print = 1:5,
  nfeatures.print = 30,
  reduction.key = "PC_",
  seed.use = 42,
  ...
) {
  # 这个函数干啥的？见后文
  # 做 features 验证，最后返回 scale.data 的子集
  data.use <- PrepDR(
    object = object,
    features = features,
    verbose = verbose
  )

  # 有传递给针对 default 的同名函数
  reduction.data <- RunPCA(
    object = data.use,
    assay = assay,
    npcs = npcs,
    rev.pca = rev.pca,
    weight.by.var = weight.by.var,
    verbose = verbose,
    ndims.print = ndims.print,
    nfeatures.print = nfeatures.print,
    reduction.key = reduction.key,
    seed.use = seed.use,
    ...

  )
  return(reduction.data)
}











()
./seurat-4.1.0/R/dimensional_reduction.R:2271:PrepDR <- function(
准备降维：常规检查步骤。


# Prep data for dimensional reduction
#
# Common checks and preparatory steps before running certain dimensional
# reduction techniques
#
# @param object        Assay object
# @param features  Features to use as input for the dimensional reduction technique.
#                      Default is variable features
# @ param verbose   Print messages and warnings
#
#
PrepDR <- function(
  object,
  features = NULL,
  slot = 'scale.data',
  verbose = TRUE
) {
  #(A1) 如果HVG长度为0，且 没有提供 features 参数，则报错
  if (length(x = VariableFeatures(object = object)) == 0 && is.null(x = features)) {
    stop("Variable features haven't been set. Run FindVariableFeatures() or provide a vector of feature names.")
  }
  #(A2) 获取 assay@scale.data
  data.use <- GetAssayData(object = object, slot = slot)
  # 如果行数为0 且 是 scale.data 则报错
  if (nrow(x = data.use ) == 0 && slot == "scale.data") {
    stop("Data has not been scaled. Please run ScaleData and retry")
  }

  #(A3) features 如果为空，则默认值为 HVG
  features <- features %||% VariableFeatures(object = object)
  # 保留在 scale.data 中的基因名，取唯一值
  features.keep <- unique(x = features[features %in% rownames(x = data.use)])
  # 如果保留的长度 小于 输入参数
  if (length(x = features.keep) < length(x = features)) {
  	# 找出差集：不在keep中的，给出警告
    features.exclude <- setdiff(x = features, y = features.keep)
    if (verbose) {
      warning(paste0("The following ", length(x = features.exclude), " features requested have not been scaled (running reduction without them): ", paste0(features.exclude, collapse = ", ")))
    }
  }
  # 覆盖输入参数
  features <- features.keep


  #(A4) 求 sd: 如果是同一个值，则sd==0，弃掉
  if (inherits(x = data.use, what = 'dgCMatrix')) {
    features.var <- RowVarSparse(mat = data.use[features, ])
  }
  else {
    features.var <- RowVar(x = data.use[features, ])
  }
  features.keep <- features[features.var > 0]
  # 如果长度小于 更新后的输入参数，则找到var=0的基因并警告
  if (length(x = features.keep) < length(x = features)) {
    features.exclude <- setdiff(x = features, y = features.keep)
    if (verbose) {
      warning(paste0("The following ", length(x = features.exclude), " features requested have zero variance (running reduction without them): ", paste0(features.exclude, collapse = ", ")))
    }
  }
  # 再次更新 输入参数
  features <- features.keep

  #(A5) 去掉na值 ?? 为什么会有NA值?
  features <- features[!is.na(x = features)]

  # 再取 scale.data 的子集
  data.use <- data.use[features, ]

  return(data.use)
}










() 最难懂的一个函数，需要奇异值分解的背景知识
./seurat-4.1.0/R/dimensional_reduction.R:847:RunPCA.default <- function(


#' @param assay Name of Assay PCA is being run on
#' @param npcs Total Number of PCs to compute and store (50 by default)
#' @param rev.pca By default computes the PCA on the cell x gene matrix. Setting
#' to true will compute it on gene x cell matrix.
#' @param weight.by.var Weight the cell embeddings by the variance of each PC
#' (weights the gene loadings if rev.pca is TRUE)
#' @param verbose Print the top genes associated with high/low loadings for
#' the PCs
#' @param ndims.print PCs to print genes for
#' @param nfeatures.print Number of genes to print for each PC
#' @param reduction.key dimensional reduction key, specifies the string before
#' the number for the dimension names. PC by default # 类似 PC_1 的前缀
#' @param seed.use Set a random seed. By default, sets the seed to 42. Setting
#' NULL will not set a seed.
#' @param approx Use truncated singular value decomposition to approximate PCA #奇异值分解 近似 PCA
#'
#' @importFrom irlba irlba
#' @importFrom stats prcomp #主力函数
#' @importFrom utils capture.output
#'
#' @rdname RunPCA
#' @concept dimensional_reduction
#' @export
#'
RunPCA.default <- function(
  object,
  assay = NULL,
  npcs = 50,
  rev.pca = FALSE,
  weight.by.var = TRUE,
  verbose = TRUE,
  ndims.print = 1:5,
  nfeatures.print = 30,
  reduction.key = "PC_",
  seed.use = 42,
  approx = TRUE,
  ...
) {
  #(A1) 如果seed.use非空，则设置随机数种子
  if (!is.null(x = seed.use)) {
    set.seed(seed = seed.use)
  }
  
  #(A2) 【跳过】如果 rev.pca ==T，则 在 gene x cell 上计算PCA。非默认过程。
  if (rev.pca) {
  	# 最小值: npcs, 列数-1
    npcs <- min(npcs, ncol(x = object) - 1)
    # 奇异值分解
    pca.results <- irlba(A = object, nv = npcs, ...)
    # 变异的总和
    total.variance <- sum(RowVar(x = t(x = object)))
    # 特征根 / 根号(行数-1)
    sdev <- pca.results$d/sqrt(max(1, nrow(x = object) - 1))

    # 左矩阵u是 gene 的loading
    # 如果 使用 var 权重
    if (weight.by.var) {
      # 则左矩阵u * 特征值对角矩阵
      feature.loadings <- pca.results$u %*% diag(pca.results$d)
    } else{
      # 否则，直接使用原 矩阵
      feature.loadings <- pca.results$u
    }

    # 右矩阵v是 cell embeddings
    cell.embeddings <- pca.results$v
  }

  #(A3) 【默认】如果 rev.pca ==F，则 在 cell x gene 上计算PCA。
  # A=u.B.v^T
  else {
    total.variance <- sum(RowVar(x = object))
    # (B1) 【默认】如果近似求解
    if (approx) {
      npcs <- min(npcs, nrow(x = object) - 1) # 至少是 行数-1
      pca.results <- irlba(A = t(x = object), nv = npcs, ...) #转置后，奇异值分解
      feature.loadings <- pca.results$v #右矩阵v是 gene 的loading
      sdev <- pca.results$d/sqrt(max(1, ncol(object) - 1)) #特征根 / 根号(列数 - 1)
      # 左奇异矩阵 u 是 cell embeddings
      if (weight.by.var) {
      	# 如果需要 var 权重
        cell.embeddings <- pca.results$u %*% diag(pca.results$d)
      } else {
        cell.embeddings <- pca.results$u
      }
    # (B2) 【跳过】如果 精确求解
    } else {
      npcs <- min(npcs, nrow(x = object)) # 至少是行数
      pca.results <- prcomp(x = t(object), rank. = npcs, ...) #转置后，做PCA
      # gene loadings
      feature.loadings <- pca.results$rotation
      # 每个PC的解释的误差
      sdev <- pca.results$sdev
      # 如果需要 var 权重，则直接使用 x
      if (weight.by.var) {
        cell.embeddings <- pca.results$x 
      } else {
      	# 如果不需要 var 权重，则去掉它们
      	# 什么情况下去掉 var 权重?为什么这样去? //todo
        cell.embeddings <- pca.results$x / (pca.results$sdev[1:npcs] * sqrt(x = ncol(x = object) - 1))
      }
    }

  }
  # 小结：以上默认走的是 A3-B1。
  # 使用近似求解，部分特征值，使用函数 irlba 求SVD 分解。
  # 使用 右奇异矩阵v 作为旋转矩阵，使用 左矩阵u . 特征根对角矩阵 作为有var权重的 cell embeding 也就是点的坐标。
  
  # 基因转换矩阵的行列
  rownames(x = feature.loadings) <- rownames(x = object)
  colnames(x = feature.loadings) <- paste0(reduction.key, 1:npcs)
  # 细胞位置的坐标
  rownames(x = cell.embeddings) <- colnames(x = object)
  colnames(x = cell.embeddings) <- colnames(x = feature.loadings)

  #创建对象
  reduction.data <- CreateDimReducObject(
    embeddings = cell.embeddings,
    loadings = feature.loadings,
    assay = assay,
    stdev = sdev,
    key = reduction.key,
    misc = list(total.variance = total.variance)
  )

  # 如果要输出信息
  if (verbose) {
  	# 这是一个神奇的函数，能获取其他函数的输出
    msg <- capture.output(print(
      x = reduction.data,
      dims = ndims.print,
      nfeatures = nfeatures.print
    ))
    # 对该输出
    message(paste(msg, collapse = '\n'))
  }
  # 返回降维对象
  return(reduction.data)
}


